clear all;
format long g;


 %% PRE-SET parameters
r      = log(1.05)/12;  
alpha  = .64;
lambda = .25;       % steady state target
En_firm=  20;       % avg firm size


 % Calibration %%%%%
printinactionrate = 27;  
printcolumn = (printinactionrate==10)*3 + (printinactionrate==17)*4 + ...
              (printinactionrate==27)*6 + (printinactionrate==38)*7;
% We will use 'printcolumn' to print results below.
%
parameter_vector = load(['/home/parameters_' num2str(printinactionrate) '/param_end_replication.mat']);
beta = parameter_vector.param_end.beta;
C = parameter_vector.param_end.K;
c = parameter_vector.param_end.c;
s = parameter_vector.param_end.s;
sigma = parameter_vector.param_end.sigma;


 


 %% Solve steady state -> distribution(m) and policy fncs
solve_main;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
 %% Simulation of individual employers
 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cores = 16;
parpool(cores,'IdleTimeout',Inf); 
 %
 %
 %% Final decisions before simulation: 
 % A) Set max growth in Net Expansion region.
Utol     = 1.75e-4;	    % This truncates super convex portion of hiring schedule.
 %
 % B) Sample size and time step:
rng(502);
NumFirms = cores*floor(5000000/cores);  %      
numfirms = NumFirms/cores;
daysteps = 12*24;   % z*24 means z steps per hour, or time step is 60/z minutes.
mosteps  = daysteps*26;
dt       = 1/mosteps;
months   = 2*12;	
NumPrds  = mosteps*months; 
 %
 % C) Set parameters of law of motion of n* = (alpha*x)^(1/(1-alpha)).
mustar = (1/(1-alpha))*mulog*dt;
sdstar = (1/(1-alpha))*sqrt(dt)*sigma;
 %
 %
 %% Initialize variables that will be saved.
 %  1) Age, Employment, Output, and Wages
InitialEmp    = zeros(numfirms, cores);             % Save initial (start of simulation) employment.
EmpEndofmo    = zeros(numfirms, cores, months );    % End-of-MONTH employment (final time step of month).
WageEndofmo   = zeros(numfirms, cores, months );    % End-of-Month Wage rate (final time step of month).
Output1mo     = zeros(numfirms, cores, months );    % Cumulates output through month.
 %
 % 2) Monthly net growth, flows, and vacancies
EmpGrowth1mo  = zeros(numfirms, cores, months );    % DFH net growth
LayoffRate1mo = zeros(numfirms, cores, months );    % DFH Rates (total monthly flow / avg employment across month)
HireRate1mo   = LayoffRate1mo;                      % ...
QuitRate1mo   = LayoffRate1mo;
HireTot1mo    = zeros(numfirms, cores, months );    % This is cumulated hires through month.
QuitTot1mo    = HireTot1mo;                         % cumulated quits 
LayoffTot1mo  = HireTot1mo;                         % cumulated layoffs
EEHireTot1mo  = HireTot1mo;                         % cumulated EE hires (quits *to* the firm) through month.
 %
 %
 %
 %% To initialize firm size, will use map from m to E[n | m]
unifm   = linspace(mL,mU,1000)';
Gf_fit_all = interp1(grid_all,Gf_all,unifm,'linear','extrap');
G_fit_all  = interp1(grid_all,G_all, unifm,'linear','extrap');
dGf_all = diff(Gf_fit_all);
dG_all  = diff(G_fit_all);
nsize0  = En_firm*dG_all./dGf_all;
mgrid0  = .5*unifm(1:end-1) + .5*unifm(2:end); 
 % Below, interpolate fnc. nsize0(mgrid0) to recover initial ns given initial ms
 %
 %
parfor i = 1:cores
    format long g;
    
     % Initial marg. product, firm size, and frictionless size (n*)
    mlag  = min( grid_all(end),  max(grid_all(1),  interp1(Gf_all, grid_all, rand(numfirms,1),'linear','extrap') )); 
    nsim  = interp1(mgrid0, nsize0, mlag,'linear','extrap');  
    nsim  = En_firm*nsim/mean(nsim); 
    %mlag  = min( grid_all(end),  max(grid_all(1),  interp1(G_all, grid_all, rand(numfirms,1),'linear','extrap') )); 
    %nsim  = En_firm*ones(numfirms,1);
    InitialEmp(:,i) = nsim;
    lastmo_nsim = nsim;
    nstarlag = (mlag.^(1/(1-alpha))) .* nsim;
     %    
     % 
     % Initialize output, wage, flows.
    rOutput1mo    = zeros(numfirms,1);
    rHireTot1mo   = zeros(numfirms,1);
    rQuitTot1mo   = rHireTot1mo;
    rLayoffTot1mo  = rHireTot1mo;
    rEEHireTot1mo  = rHireTot1mo;
     %
     %
     %
    for mo = 1:months
        %tic
        for t = 1:mosteps
             %
             %% Draw TFP innovation:
            dz    = randn(numfirms,1);
            dlognstar = mustar + sdstar*dz;
            nstar = nstarlag.*exp(dlognstar);
            xt    = (1/alpha)*(nstar.^(1-alpha));	% n* = (alpha*x)^(1/(1-alpha))
             %
             %
             %% OPTIMAL POLICY
            nt = zeros(numfirms,1);
            mt = zeros(numfirms,1);
            growthn  = zeros(numfirms,1);
            layoffrate = zeros(numfirms,1);
            hirerate = zeros(numfirms,1);
            quitrate = zeros(numfirms,1);
            seprate  = zeros(numfirms,1);
             %    
             %
             % I) Layoffs:
             %%%%%%%%%%%%%%%%%%%%%%%%%%%
             % After (all) quits, define 
            nhat  = nsim - s*lambda*nsim*dt;
             %
             % What if didn't hire? Marginal product would be
            mhat  = alpha*xt.*(nhat.^(alpha-1)); 
            fire1 = (mhat <= mL);       % Lay off if mpl still < mL *after* all quits.
            fire  = find(fire1 ==1);
             %
             % Net growth: Lower n s.t. mpl = mL -->
            mt(fire) = mL; 
            nt(fire) = (alpha*xt(fire)./mL).^(1/(1-alpha));
            growthn(fire) = (nt(fire) - nsim(fire))./nsim(fire);     % recall, nsim = lagged n.
             %
             % Turnover. Note layoff rate defined by layoffs == nhat - nt == nsim*LAYOFF RATE*dt.
            layoffrate(fire) = (nhat(fire) - nt(fire))./nsim(fire) / dt;
            quitrate(fire) = s*lambda;
            seprate(fire) = s*lambda + layoffrate(fire); 
             %
             %
             % II) Natural wastage
             %%%%%%%%%%%%%%%%%%%%%
            nw1 = (mL < mhat) .* (mhat <=mM);
            nw  = find( nw1 ==1);
            nt(nw) = nhat(nw);
            mt(nw) = mhat(nw);
             % Net growth:
            growthn(nw)  = -s*lambda*dt;
             % Turnover:
            quitrate(nw) = s*lambda;
            seprate(nw) = s*lambda;  
             %
             %
             % III) Full Replacement:
             %%%%%%%%%%%%%%%%%%%%%%%%
             % Net growth = 0 --> marg. product is 
            m_fr  = alpha*xt.*(nsim.^(alpha-1)); 
             %
             % If mM < m_fr < mR, then unambiguously in FR region.
            fr_a  = (mM < m_fr).*(m_fr <=mR); 
             %
             % However, w/ discrete time steps, it may happen that 
             % m_fr <mM < mhat. This says that mpl if fully replace  
             % (and thus, n(t) = lagged n) is in the full attrition region 
             % after x arrives (m_fr <mM). BUT if allow full attrition, mpl  
             % drifts above mM threshold (mM< mhat). I assume these employers
             % fully replace. This case accounts for very small share of 
             % employers -- see below.
            fr_b  = (m_fr < mM).*(mM < mhat);
            fr1   = fr_a  +  fr_b;
            fr    = find(fr1 ==1);
            mt(fr)= m_fr(fr);
             % Size 
            nt(fr)= nsim(fr);
             %
             % Turnover:
            quitrate(fr) = interp1(grid_fr, quit_fr, m_fr(fr), 'linear', quit_fr(1));
             % Since a handful of m_fr < mM==grid_fr(1), specify that quit
             % rate = quit_fr(1) in these cases.
            seprate(fr)  = quitrate(fr);
            hirerate(fr) = seprate(fr);
             %
             %
             % IV) Net Expansion:
             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             % MPL before any emp. adjustment:
            ne1   = (mR <m_fr);     % **if** you were to replace and still mpl is too big, then hire.
            nepop = sum( ne1 );
            if nepop >0
                ne    = find(ne1 ==1);
                 % First, interpolate quit and hiring rates @ mpl BEFORE adjustment. 
                 % Then compute implied firm size and mpl. 
                mhat_ne     = alpha*xt(ne).*(nsim(ne).^(alpha-1));  
                quitrate(ne)= interp1(grid_ne, quit_ne,  max(grid_ne(1), min(mhat_ne, grid_ne(end))), 'linear');
                hirerate(ne)= interp1(grid_ne, eta_ne,   max(grid_ne(1), min(mhat_ne, grid_ne(end) - Utol)), 'pchip');
                 % Subtract Utol off grid_ne(end) to avoid the steepest
                 % portion of hiring rate fnc. (where numerical error 
                 % most likely).
                seprate(ne)  = quitrate(ne);
                growthn(ne)  = ( hirerate(ne) - seprate(ne) )*dt;
                nt(ne)       = nsim(ne)  + growthn(ne).*nsim(ne);
                mt(ne)   = alpha*xt(ne).*(nt(ne).^(alpha-1));    % mpl after adjustment
                Avgm_ne  = mean(mt(ne));                         % avg. mpl among net expanders.
                 %                    
                 % Check that interpolation executed correctly.
                quitnans = sum(isnan(quitrate(ne))) + sum((quitrate(ne) <0));
                hirenans = sum(isnan(hirerate(ne))) + sum((hirerate(ne) <0));
                if quitnans >0 || hirenans >0
                    disp('Error in calculating flows in NE region');
                    %pause;
                end                 
            end
             %
             %
             %% Check for errors in calculation of policy 
             %  Check #1:
            regimes = fire1 + nw1 + fr1 + ne1;
            if sum(regimes) ~= numfirms
                disp('Error: Multiple regimes or failure to assign any regime');
                %pause;
            end
             % Check #2:
            sumfr1   = sum(fr1);
            sumfr_b  = sum(fr_b);
            shrfr_b  = sumfr_b / sumfr1
            if shrfr_b > 0.001
                disp('Error: Failure in assigning attrition v. replacement regime');
            end
             %  Check #3:
            indgrow_fire = (growthn(fire) >0);
            if nepop >0
                indgrow_ne   = (growthn(ne)   <0);
            else 
                indgrow_ne   = 0;
            end
            if sum(indgrow_fire) >0 || sum(indgrow_ne) >0
                disp('Error: Net growth is inconsistent with assigned regime');
                %pause;
            end
             %
             %
             %% EE share of hires (these are quits *to* the firm) 
            [~, pos] = sort(mt);    
            simG     = zeros( numfirms, 1 );
            simG(pos)= cumsum(nt(pos)) / sum(nt);       % Emp-weighted c.d.f. of marg. product
            uratesim = ( L - mean(nt) ) / L;            % ( L produced in solve_moments.m )
            eeshrhires = s*(1-uratesim)*simG ./ (uratesim + s*(1-uratesim)*simG);
            eehirerate = eeshrhires.*hirerate; 
             %
             %
             %% Cumulate/update major aggregates. 
            rHireTot1mo   = rHireTot1mo    + hirerate.*nsim*dt;              
            rQuitTot1mo   = rQuitTot1mo    + quitrate.*nsim*dt;              
            rEEHireTot1mo = rEEHireTot1mo  + eehirerate.*nsim*dt;                  
            rLayoffTot1mo = rLayoffTot1mo  + layoffrate.*nsim*dt;            
             %
             %  Next up - Age, avg. size, monthly output and wages: 
            apl        = mt / alpha; 
            wage       = (beta/(1-beta*(1-alpha)))*mt + omega0;
            rOutput1mo = rOutput1mo + dt*apl.*nt;
             %
             %
             %% Reset lagged, n, n*, and m:
            nsim = nt;
            mlag = mt;
            nstarlag = nstar;
        end     %%%  END time steps WITHIN MONTH %%%
         %
         %
         % Report progress:
        disp( [i, mo] );  
         %    
         %
         % Size, output and wages 
        EmpEndofmo(:,i,mo)  = nsim;     
        WageEndofmo(:,i,mo) = wage;
        Output1mo( :,i,mo)  = rOutput1mo; 
        rOutput1mo   = zeros(numfirms,1);   % Reset    
         %
         %
         % Monthly flows
        HireTot1mo(:,i,mo)     = rHireTot1mo;
        QuitTot1mo(:,i,mo)     = rQuitTot1mo;
        LayoffTot1mo(:,i,mo)   = rLayoffTot1mo;
        EEHireTot1mo(:,i,mo)   = rEEHireTot1mo;
        EmpGrowth1mo( :,i,mo)  = (nsim - lastmo_nsim) ./ (0.5*(nsim + lastmo_nsim));
        LayoffRate1mo(:,i,mo)  = rLayoffTot1mo  ./ (0.5*(nsim + lastmo_nsim));
        QuitRate1mo(:,i,mo)    = rQuitTot1mo  ./ (0.5*(nsim + lastmo_nsim));
        HireRate1mo(:,i,mo)    = rHireTot1mo  ./ (0.5*(nsim + lastmo_nsim));
         %
         % Reset:
        rHireTot1mo    = zeros(numfirms,1);
        rQuitTot1mo    = rHireTot1mo;
        rLayoffTot1mo  = rHireTot1mo;
        rEEHireTot1mo  = rHireTot1mo;
        lastmo_nsim    = nsim;
        
    end     %%%  END MONTH %%%
end     %%%  END CORES %%%
 %
 %
 % Reshape to compute moments.
Emp0           = reshape(InitialEmp, NumFirms, 1);
Output1mo      = reshape(Output1mo,      NumFirms, months);
EmpEndofmo     = reshape(EmpEndofmo,     NumFirms, months);
WageEndofmo    = reshape(WageEndofmo,    NumFirms, months);
%
EmpGrowth1mo   = reshape(EmpGrowth1mo,   NumFirms, months);
LayoffRate1mo  = reshape(LayoffRate1mo,  NumFirms, months);
HireRate1mo    = reshape(HireRate1mo,    NumFirms, months);
QuitRate1mo    = reshape(QuitRate1mo,    NumFirms, months);
LayoffTot1mo   = reshape(LayoffTot1mo,   NumFirms, months);
EEHireTot1mo   = reshape(EEHireTot1mo,   NumFirms, months);
HireTot1mo     = reshape(HireTot1mo,     NumFirms, months);
QuitTot1mo     = reshape(QuitTot1mo,     NumFirms, months);


 %
 % DFH Hockey Stick
sim_hockey; 
sheetname = ['Fig 5 - ' num2str(printinactionrate)];
writematrix(hockeyMAT, '/home/Output/ExtraChains.xlsx', 'Sheet', sheetname, 'Range','A2:G62');
% hockeyMAT    = [min of bin, max of bin, midpt of bin, Hiring rate, Quit rate, Sep. rate, Layoff rate];  


 %
 % Employment Growth
 % (This pulls from variables made in sim_hockey).
EEmpgrowth    = (qtAvgEmpVec'/qtTotAvgEmp)*qtEmpGrowthVec;
VarEmpgrowth  = (qtAvgEmpVec'/qtTotAvgEmp)*((qtEmpGrowthVec - EEmpgrowth).^2);
StdvEmpgrowth = sqrt(VarEmpgrowth);
disp('Std dev. of Quarterly Employment growth =');
disp( StdvEmpgrowth );
rangename = append('B', num2str(printcolumn));
writematrix(StdvEmpgrowth, '/home/Output/Chains.xlsx', 'Sheet', 'Table 2, Table 4', 'Range', rangename);

 %
 % Misallocation (dispersion of APL)
sim_apl;
StdvlogApl = sqrt(var_logApl);
disp('Std dev of log APL');  disp( StdvlogApl );
disp('Autocorrelation of log APL');  disp(rho_logApl);
disp('Std dev. of APL growth: Inactive, Adjusters');
StdvAplgrowth = sqrt( [VInAPLGrowth12mo; VAdjAPLGrowth12mo] );
disp( StdvAplgrowth );
%
rangename = append('C', num2str(printcolumn));
writematrix(StdvlogApl, '/home/Output/Chains.xlsx', 'Sheet', 'Table 2, Table 4', 'Range',rangename);
rangename = append('D', num2str(printcolumn));
writematrix(rho_logApl, '/home/Output/Chains.xlsx', 'Sheet', 'Table 2, Table 4', 'Range',rangename);
sheetname = ['Table F2 - ' num2str(printinactionrate)];
writematrix(StdvAplgrowth, '/home/Output/ExtraChains.xlsx', 'Sheet', sheetname, 'Range','B2:B3');

 %
 % Poaching index calculations:
AvgEEShare = sum(EEHireTot1mo,2) ./ sum(HireTot1mo,2);    % NumFirms by 1
sim_poach;
sheetname = ['Fig E - ' num2str(printinactionrate)];
writematrix(PoachBinMAT, '/home/Output/ExtraChains.xlsx', 'Sheet', sheetname, 'Range','A3:D12');
% PoachBinMAT   = [percentile,  Wage rank, EU , EE ];


delete(gcp('nocreate'));
